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NONLINEAR TRANSIENT PROBLEMS USING STRUCTURE - 
COMPATIBLE HEAT TRANSFER CODE 


Gene Hou, Yang Wang and Rufeng Liu 

Department of Mechanical Engineering 
Old Dominion University 
Norfolk, VA 23454 

The report documents the recent effort to enhance a transient linear heat transfer 
code so as to solve nonlinear problems. The linear heat transfer code was originally 
developed by Dr. Kim Bey of NASA Largely and called the Structure-Compatible Heat 
Transfer (SCHT) code. The report includes four parts. The first part outlines the 
formulation of the heat transfer problem of concern. The second and the third parts give 
detailed procedures to construct the nonlinear finite element equations and the required 
Jacobian matrices for the nonlinear iterative method, Newton-Raphson method. The final 
part summarizes the results of the numerical experiments on the newly enhanced SCHT 
code. 

I. FORMULATION 

The system of the governing differential equations of a heat transfer problem in a 
homogeneous material can be expressed as 


and 
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where the temperature, T ( x , t), is the only unknown. Equation (1) represents an initial- 
boundary value problem with Eq. (2) and Eq. (3) being the temperature and the heat flux 
boundary conditions, respectively, and Eq.(4), the initial conditions. It is assumed that the 
heat source, Q, the prescribed temperature,/, and the flux, q s , and the initial value, g, are 
all with proper regularity. In case of material nonlinearity, the coefficients in Eq.(l), the 
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specific heat, pcij ) and the thermal conductivity, ( T) are assumed to be functions of 

temperature. 

The weak form of the above equations can be derived based upon the 
discontinuous Galerkin method [1] for a time interval, /„ = [ t n .\, t n ) as 
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where the testing function, v = 0 on dQ u and the temperature be a combination of the 
unknown function, u, and an arbitrary function, w. The testing function should satisfy the 
boundary condition, w = /, on dQ. u . Consequently, u = 0 on dQ u . Note that u is the 
only unknown in Eq. (5). 

The symbols of T~_ t and T n + _, are the limits of temperatures approach to the point, 
t n X , from its right and left sides. It should be noted that the continuity of the temperature, 

T, is not enforced pointwisely throughout the entire time domain, rather it is enforced at 
the ends of the time interval in a weighted manner. Furthermore, the coefficients, pc and 
Kij, also experience discontinuity at the ends of time intervals, as they are functions of the 
temperatures, which are allowed to be discontinuous at those points. 


II. FINITE ELEMENT EQUATIONS 

For simplicity, one will start the derivation with an assumption that the problem is 
homogeneous. That is,/= 0 and T = it. 

The weak equation defined in Eq. (5) can be discretized as a summation of 
integrals over individual "finite elements". In this study, the computational domain of 
such a finite element, V/, is defined as a product of a spatial domain and a time interval, 
Q n Xl n . Asa result, the discretized equations for time interval, /„, become 


2£fkfr v+ (*- v * r )- v * v 


A 

cbc dt + ^^pcT^vdx 

J Q " 


2 



( 6 ) 


= ?£, (L. Qvdx ) dt + XL. ^ r .'-' v<4c + X£, L„. v. vdsdl 


where v is the testing function. The tasks in constructing the equivalent matrix equation 
of Eq. (6) involve selection of the interpolation functions, approximation of nonlinear 
material properties and integration of elemental matrices. 

Interpolation Functions 

Following a similar derivation given in [1], a hp-wt rsion of interpolation function 
is introduced here to approximate the solution in the kth element space domain and the (n 
- l)th time interval as 


T(x,t) = X(x,i)a nk 


where X ( x , t) is the interpolation function vector and a n k is the "nodal value" vector. 

More specifically, the interpolation function can be spanned as a tensor product of 
functions with separable variables as 

X(x,t)= L(x,y)®y/(t)®<p(z) 


or in the index form. 
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The subscript, m, in the above equation follows the relation, m = (i-l)ps + (J- 
1) s + k where p and s are the ranges of subscripts, j and k. In the case of a triangular 
element, the interpolation function, L„ can be simply specified as an area coordinate, 

L { = a t +b i x + c i y , i = 1,2,3 

The interpolation function through the thickness. A, is a polynomial of a 
hierarchical basis. In this study, <p(z) is given as 


3 




1 - 2 - 


) 



1 + 2 - 


J 


' V2(2 j - 1) 


r ( 7 \ 

2 - 

v v A y 




r 7 V\ 

2 - 

v A yy 


7 = 2,..., p 


! K*i 



where z ranges from A/2 to A/2 , and Pj is a Legendre polynomial of degree j. 

The interpolation function for the time domain, y/(t), uses a power polynomial in 
time. That is, 


Vj j = 0,...,s 


Nonlinear Approximation 

It is assumed that in this study, the relations between the material properties, pc 
and Kij, and the temperature are defined in a tabulated form, presented as a result of 
experiments. Therefore, the material properties can not be explicitly specified as 
functions of position and time as required by integration. An approximation is thus 
introduced to overcome such a difficulty. 

A single hp-e lement is first divided into regular subdomains. The values of the 
material properties at the vertices of each subdomain are read out from the given material 
table based upon the values of the temperature found at those points. The values of the 
material properties at elsewhere in a subdomain are then obtained through linear 
interpolation. In this way, the material properties are explicitly approximated as functions 
of position and time throughout the problem domain. 

As an example, the spatial and the time domain of an element can be divided into 
subdomains as shown in Fig. 1. and the material property, say pc, can be interpolated 
linearly in a typical subdomain bounded by [zi , Z 2 ] x [ ti , r? ] as 
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pc(x , t) = ( l ( jc , y) ® N t (/) ® N 2 (z)Y c 

where c is the vector of the nodal values of the material properties that are found through 
the material table and the nodal values of temperature. The area coordinates are in this 
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case, a nature choice for linear interpolation functions in x-y plane, L (x, y), because of a 
triangular domain. Typical linear shape functions can be used to interpolate N, and N z as 
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As a result, the vector c has 12 components, which gives an approximation of the 
temperature field in the subdomain of concern as 


£ (^1 5 ^2 ’ ^3 ’ ^4 » ^5 » ^6 ’ *'7 ’ ^8 » ^9 ’ ^ 10 ’ ^ i 1 » ^12 ) 
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where the first subscript of the notation, cp, is associated with the area coordinates, the 
second with the vertices of a time subdomain and the third with the vertices of a z 
subdomain. Thus, Cy*, is read from the given material table based upon the value of the 
temperature evaluated at area nodal point i, time vertex j and z vertex k. As an example, 
the material property, c/ 22 » is obtained at the point ( xi, yi, Zi ) x t 2 as 
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Matrix Evaluation 


Once the interpolation functions are selected and the nonlinear material properties 
are approximated as explicit functions of position and time, one can proceed to integrate 
the terms in Eq. (6) to construct the equivalent matrices. The resultant finite element 
matrix equation for the time interval, I n , can then be expressed as 


(C„ + K, + M, Xfl, } = {q , } + M.„ {a„ J+ {*„ } (7) 

where the subscript, n, indicates that the associated matrix or vector is evaluated with 
functions defined in time interval, /„ . Note that each term in Eq. (7) is corresponding to 
an integral in Eq. (6). As an example, the capacitance matrix is given as 


- 
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For a specific element k which includes a triangular slab in space and /„ in time, 
the elemental capacitance matrix with constant pc can be further spelled out as 

C . = f" If pcX^—dx\t 

"•* dt J 

= t, l y) ® v{* ) ® <p(z)\ L (x, t) ® x//(t) ® ^(2)) dxdt 

J / T A 

= f LlJ dxdy® P" \l / — dt® \\ pc<p(p T dz 
JA j '"-' dt J -i 

where A is the area of the triangular element. In the case of nonlinear pc , the elemental 
C matrix is a function of temperature, T, and can be written as 
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Thus, the elemental capacitance matrix is a summation of the subdomain capacitance 
matrices. The more subdomains divided in each element, the more computational efforts 
are required. 
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III. SOLUTION PROCEDURES 


The matrix equation derived above is based upon the discontinuous Galerkin’s 
method and will be solved in a time-marching fashion. In other words, the matrix 
equation of Eq. (7) will be solved for a time interval at a time, with a n as unknown and a n . 
i as known quantities. The value, ao , is given as the initial condition. The nonlinear 
version of Eq.(7) is rewritten below for future reference 

[C n k ) + K n {a J + M n (a n )\a n }= {q n }+ [M„_, (fl H )\a n _ l }+ {b n } 


The above equation can be solved by either the fixed point iteration or the Newton- 
Raphson’s method. 

The recursive formula for the fixed point iteration is simply given as 

[C, (<-')+ K, (<")+ M. (<" )k}= {<!. 1+ [M,., )k-, }+{*,} 


where the superscript, i , denotes the iteration number. The recursive formula for the 
Newton-Raphson’s method, however, is more complicated and involves Jacobian 
matrices, J c , J k and J m as 
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and the solution is updated by 


a -a 


i " 1 +A al 


(10) 


In the above equation, A a' n is the improvement of the solution and R' n 1 is the residual of 
the nonlinear equation at the i-l iteration which is defined as 


{/?;' } = [C n [a ) + K n (a '' 1 ) + M n (a l ' 1 }- {q n }[ M n _ x (a„_, )]{«„.,}- {b n } 


Moreover, the Jacobian matrices, J c , 7* and 7m are obtained by differentiating the 
coefficient matrices C„ , K n and M n with respect to the unknown vector a n , respectively. 
This is usually accomplished at the element level. As an example, the Jacobian matrix of 
the element capacitance matrix is obtained for a typical element ( n, k ) as 
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where dC n k / dT t is obtained by using Eq. (8) as 
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where 3c iyit / 97 yt is the derivative of pc with respect to the temperature at * , y t , Zj , and 
t k . This derivative is based upon the relation between the material properties and the 
temperature described in the material table. Since the temperature T jjk is given as 

T, Jk = [£(*, , y, ) ® ¥(*, ) ® <p(z k ^ T , 


the derivative, dT jjk / dT t , is the coefficient of T t in the above equation. 

In the Newton-Raphson’s iteration, Eq. (9) is first solved for the solution 
improvement A a n _k which is then updated based upon Eq. (10) to get a better solution. A 
Newton-Raphson’s iteration will be restarted with the new solution until the error residual 
is reduced to an accepted level. 


IV. NUMERICAL RESULTS 

Several examples have been studied to verify the accuracy of the solution 
procedure, designed to solve the nonlinear finite element equations derived by the 
discontinous Galerkin’s method for the heat transfer problems. To start the example 
problems, an exact solution of temperature is first selected, which gives zero temperature 
on the entire surface of the domain. That is, the example problems have zero boundary 
temperature. The heat source term, Q, that generates the prescribed temperature, can be 
obtained by substituting the given temperature function into the governing differential 
equation of Eq. (1). The material properties are assumed to be linear in temperature and 
in the form of (a T + b ). 

The convergence criterion used in this study is defined that every component of 
the unknown vector, a, should satisfy the condition 
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IV. 1 One Dimensional Problem 

The exact temperature field is given as 

T(x,t) = t^n-<p n (rj) 

rt =2 

The fifth order polynomial is used in the finite element equation to approximate 
the solution through thickness. The entire thickness is divided into three subdomains, in 
each of which the material properties will be expressed as a linear function of position 
and time. 

Case 1: Effect of Nonlinearity 

The constants in the assumed relation between the material properties and the 
temperature will be varied to study the effect of nonlinearity on solution accuracy. In 
this study, the material properties are assumed to be in the form of (T+ b). The values 
of b ranges from 3 to 15. Since the maximal value of the temperature in this problem is 
around 2.5, the problem will experience less nonlinearity with a higher value of b. 

Figure 2 shows the temperature solutions with various values of b. It is shown 
that, under the current discretization, the solution procedure is not able to reach at the 
correct solution when the nonlinearity is high. The results of numerical experiments also 
show that using the converged solution of the problem with less nonlinearity as the initial 
guess does not improve the solution accuracy of the problem with higher nonlinearity. 

Case 2: Approximation of Time Axis 

This case will study the effectiveness of the approximation along the time axis. 
The exact temperature field is given as a quadratic function of time as 

T(x,t)=e±„t(v) 

n=2 


Figures 3 and 4 compare the calculated temperatures with the exact solution at 
time =1 and 2 seconds. The symbols ipt, dt and nt in the figures are referred to the order 
of polynomial of time, the time interval for each element and the number of subdomains 
along the time axis, respectively. It is shown that increasing ipt and nt, and reducing dt 
do not have significant impact on the accuracy of the solution. 
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IV. 2 Two Dimensional Examples 


The domain of the problem is a 1 x 1 x 1 slab. The exact solution is given as 

T(x,t) = txy(x - l)(j - l)X n<p n (rj) 

n= 2 


The central plane of the slab is discretized into three models with 8, 32 and 128 
triangular elements, respectively. These models are shown in Figs. 5 to 7. For most of the 
cases studied, a 5th order polynomial in z and a linear polynomial in time are used for 
interpolation. The material properties are in the form of ( T + 1). The convergence 
criterion of Eq. (11) is set to be 0.001. Finally, both of the time interval and the thickness 
are discretized into three subdomains. 


Case 1: Better Approximation in Material Nonlinearity 

Figures 8 and 9 shows the comparisons of results of the 8-element model reported 
at the center of the domain for three time instances. The results in Fig.8 are obtained 
with both the z and the time intervals being discretized into 3 subdomains, while those in 
Fig. 9 discretized into 6 subdomains. The results reveal that increasing the subdomains to 
better the approximation of the material properties does not uniformly improve the 
accuracy through the time and the spatial domains. 

Similar studies are done with 32-element and 128-element models. The results are 
summarized in Figs. 10 and 11 and Figs. 12 and 13. Note that only the number of 
subdomains along the z interval is increased in these cases. Again, it is observed that 
increasing the number of points to better approximate the material nonlinearity may not 
necessarily result in an improvement in the solution accuracy everywhere. 

Comparison of Figs. 8 and 9 to Figs. 10 and 11 and to Figs. 12 and 13 reveals 
another concern that reducing the size of the h elements does not improve the solution. 
Similar study, done on cases with linear material properties, draws a similar conclusion. 
The results are shown in Figs. 14 to 16. This observation leads to the next case of study. 


Case 2: Reduction of h-size 


The temperatures at the off-centered points, 59 and 61 on Fig. 7, are calculated 
and reported in Figs. 17 and 18, respectively. It is clear that the reduction of element size 
does not warrant an improved solution uniformly throughout the time and the spatial 
domains. 
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V. SUGGESTION FOR THE FUTURE STUDY 


The initial study shows that more works need to be done in order to make the 
discontinous Galerkin’s method an effective one for nonlinear transient problems. On 
one hand, one needs to develop a theoretical base to control the numerical errors in terms 
of orders of the polynomial and size of the elements. On the other hand, one needs to 
develop an effective method to approximate the material nonlinearly and to reduce 
computational cumbersome of the Jacobean matrices. 
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Figure 2. Accuracy of Solutions with Different Degrees of Nonlinearity 
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Figure 3. Different Approximation along Time Axis; Solutions at Second 1 
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Figure 4. Different Approximation along Time Axis; Solutions at Second 2 
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Figure 8. Solutions for 8-Element Model with 3 Subdomains in z and t 
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Figure 9. Solutions for 8-Element Model with 6 Subdomains in z and t 
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Figure 10. Solutions for 32-Element Model with 3 Subdomains in z and t 
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Figure 11. Solutions for 32-Element Model with 6 Subdomains in z and 3 Subdomains in t 
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Figure 12. Solutions for 128-Element Model with 3 Subdomains in z and t 
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Figure 13. Solutions for 128-Element Model with 6 Subdomains in z and 3 Subdomains in t 
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Figure 16. Linear Transient Solutions with 128-Element Model 
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